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Abstract. Complex networks have become increasingly popular for mod¬ 
eling various real-world phenomena. Realistic generative network models 
are important in this context as they avoid privacy concerns of real 
data and simplify complex network research regarding data sharing, re¬ 
producibility, and scalability studies. Random hyperbolic graphs are a 
well-analyzed family of geometric graphs. Previous work provided empir¬ 
ical and theoretical evidence that this generative graph model creates 
networks with non-vanishing clustering and other realistic features. How¬ 
ever, the investigated networks in previous applied work were small, 
possibly due to the quadratic running time of a previous generator. 

In this work we provide the first generation algorithm for these networks 
with subquadratic running time. We prove a time complexity of -I- 

m) log n) with high probability for the generation process. This running 
time is confirmed by experimental data with our implementation. The 
acceleration stems primarily from the reduction of pairwise distance 
computations through a polar quadtree, which we adapt to hyperbolic 
space for this purpose. In practice we improve the running time of a 
previous implementation by at least two orders of magnitude this way. 
Networks with billions of edges can now be generated in a few minutes. 
Finally, we evaluate the largest networks of this model published so far. 
Our empirical analysis shows that important features are retained over 
different graph densities and degree distributions. 

Keywords: complex networks, hyperbolic geometry, generative graph 
model, polar quadtree, network analysis 


1 Introduction 

The algorithmic analysis of complex networks is a highly active research area 
since complex networks are increasingly used to represent phenomena as varied 
as the WWW, social relations, protein interactions, and brain topology [El- 
Complex networks have several non-trivial topological features: They are usually 
scale-free, which refers to the presence of a few high-degree vertices (hubs) among 
many low-degree vertices. A heavy-tail degree distribution that occurs frequently 

* This work is partially supported by SPP 1736 Algorithms for Big Data of the German 
Research Foundation (DFG) and by the Ministry of Science, Research and the Arts 
Baden-Wiirttemberg (MWK) via project Parallel Analysis of Dynamic Networks. 



in practice follows a power law |16[ Chap. 8.4], i. e. the number of vertices with 
degree k is proportional to k~'^, for a fixed exponent 7 > 0. Moreover, complex 
networks often have the small-world property, i. e. the typical distance between 
two vertices is surprisingly small, regardless of network size and growth. 

Generative network models play a central role in many complex network 
studies for several reasons: Real data often contains confidential information; it 
is then desirable to work on similar synthetic networks instead. Quick testing of 
algorithms requires small test cases, while benchmarks and scalability studies 
need bigger graphs. Graph generators can provide data at different user-defined 
scales for this purpose. Also, transmitting and storing a generative model and its 
parameters is much easier than doing the same with a gigabyte-sized network. 
A central goal for generative models is to produce networks which replicate 
relevant structural features of real-world networks [ 3 . Finally, generative models 
are an important theoretical part of network science, as they can improve our 
understanding of network formation. The most widely used graph-based system 
benchmark in high-performance computing, Graph500 [6] , is based on R-MAT m- 
This model is efficiently computable, but has important drawbacks concerning 
realism and preservation of properties over different graph sizes |12| . 

Random hyperbolic graphs, first presented by Krioukov et al. |13j , are a very 
promising graph family in this context: They yield a provably high clustering 
coefficient m, small diameter [5] and a power-law degree distribution with 
adjustable exponent. They are based on hyperbolic geometry, which has negative 
curvature and is the basis for one of the three isotropic spaces. (The other two are 
Euclidean (flat) and spherical geometry (positive curvature).) In the generative 
model, vertices are distributed randomly on a hyperbolic disk of radius R and 
edges are inserted for every vertex pair whose distance is below i?Q This family 
of graphs has been analyzed well theoretically mm and Krioukov et al. m 
show that complex networks have a natural embedding in hyperbolic geometry. 

Calculating all pairwise distances in their generation process has quadratic 
time complexity. This impedes the creation of massive networks and is likely the 
reason previously published networks based on hyperbolic geometry have been 
in the range of at most 10® vertices. A faster generator is necessary to use this 
promising model for networks of interesting scales. Also, more detailed parameter 
studies and network analytic comparisons are necessary in practice. 

Outline and contribution. We develop, analyze, and implement a fast, sub¬ 
quadratic generation algorithm for random hyperbolic graphs. 

To lay the foundation. Section introduces fundamentals of hyperbolic 
geometry. The main technical part starts with Section]^ in which we use the 
Poincare disk model to relate hyperbolic to Euclidean geometry. This allows 
the use of a new spatial data structure, namely a polar quadtree adapted to 
hyperbolic space, to reduce both asymptotic complexity and running time of the 
generation. We further prove the time complexity of our generation process to 

^ We consider the name “hyperbolic unit-disk graphs” as more precise, but we use 
“random hyperbolic graphs” to be consistent with the literature. 
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be 0((n^/^ + m) log n) with high probability (whp, i. e. > 1 — 1/n) for a graph 
with n vertices, m edges, and sufficiently large n. 

In Section]^ we add to previous studies a comprehensive network analytic 
evaluation of random hyperbolic graphs. This evaluation shows many realistic 
features over a wide parameter range. The experimental results also confirm the 
theoretical expected running time. A graph with 10^ vertices and 10® edges can 
be generated with our shared-memory parallel implementation in about 8 minutes. 
The generator will be made available in a future version of NetworKit [50] , our 
open-source framework for large-scale network analysis. Material omitted due 
to space constraints can be found in the appendix. 

2 Related Work and Preliminaries 

Related generative graph models are discussed in Section |4.3t where we also 
compare them to random hyperbolic graphs, in part by using empirical data. 


2.1 Graphs in Hyperbolic Geometry 


Kriokouv et al. m introduce the family of random hyperbolic graphs and show 
how they naturally develop a power-law degree distribution and other properties 
of complex networks. In the generative model, vertices are generated as points 
in polar coordinates {(j), r) on a disk of radius R in the hyperbolic plane with 
curvature — C^- We denote this disk with Djj. The angular coordinate 4> is drawn 
from a uniform distribution over [0, 27r]. The probability density for the radial 
coordinate r is given by [13 Eq. (17)] and controlled by a growth parameter a: 


, sinh(ar) 

“cosh(aRr- T 


( 1 ) 


For a = 1, this yields a uniform distribution on hyperbolic space within For 
lower values of a, vertices are more likely to be in the center, for higher values 
more likely at the border of Dij. 

Any two vertices u and v are connected by an edge if their hyperbolic distance 
dist^(u,r!) is below R. (Kriokouv et al. also present a more general model in 
which edges are inserted with a probability depending on hyperbolic distance. 
This extended model is not discussed here.) The neighborhood of a point (= 
vertex) thus consists of the points lying in a hyperbolic circle around it. Krioukov 
et al. show that for a/^ > i, the resulting degree distribution follows a power 
law with exponent 2 ■ a/C, + 1 [13 Eq. (29)]. Gugelmann et al. [11] analyze this 
model theoretically and prove non-vanishing clustering and a low variation of 
the clustering coefficient. Bode et al. [3 discuss the size of the giant component 
and the probability that the graph is connected [Sj- They also show [7] that the 
curvature parameter C can be fixed while retaining all degrees of freedom, we 
thus assume C = 1 from now on. The average degree fc of a random hyperbolic 
graph is controlled with the radius i?, using m Eq. (22)]. 
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(a) Native representation ([S])- 



(b) Poincare disk model. 


Fig. 1; Comparison of geometries. Neighbors of the bold blue vertex are in a 
hyperbolic respective Euclidean circle. 


An example graph with 500 vertices, R ~ 5.08 and a = 0.8 is shown in 
Figure For the purpose of illustration in the figure, we add an edge {u, v) 
where dist'H(M,u) < 0.2 • R. The neighborhood of u (the bold blue vertex) then 
consists of vertices v within a hyperbolic circle (marked in blue). 

A previous generator code with a more general edge probability and quadratic 
time complexity is available [3]. We show in Section 4.4 that for the random 
hyperbolic graphs described above, our implementation is at least two orders of 
magnitude faster in practice. 


2.2 Poincare Disk Model 

The Poincare disk model is one of several representations of hyperbolic space 
within Euclidean geometry and maps the hyperbolic plane onto the Euclidean 
unit disk Di(0). The hyperbolic distance between two points petQe G -Di(O) is 
then given by the Poincare metric |4j: 

dist«(pE, qe) = acosh (l + 2— - || 2 ^ ) ' (2) 

V (1 - lbBlr)(i - Ikisir)/ 

Figure [Tb| shows the same graph as in Figure [T^ but translated into the Poincare 
model. This model is conformal, i. e. it preserves angles. More importantly for us, 
it maps hyperbolic circles onto Euclidean circles. 


3 Fast Generation of Graphs in Hyperbolic Geometry 

We proceed by showing how to relate hyperbolic to Euclidean circles. Using this 
transformation, we are able to partition the Poincare disk with a polar quadtree 
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Algorithm 1: Graph Generation 


1 

2 

3 

4 

5 

6 

7 
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Input: n, fc, a 
Output: G = (G, E) 

R = getTargetRadius(n, fc, a); /* 

V = n vertices; 

T = empty polar quadtree of radius mapToPoincare(i?); 
for vertex v G V do 

draw (j)[v] from W[0, 2n)-, 

draw r'u[v] with density /(r) = q sinh(r)/(cosh(Qi?) — 1); 
teIv] — mapToPoincare(r-HH); 
insert v into T at {(f)[v\,rE[v])', 


Eq. (22) lEllI 


/* Eq.0 


*/ 


*/ 


9 for vertex v G V do in parallel 

10 C-H = circle around {(l>\v\,r-H[v\) with radius R\ 

11 Ce = transformCircleToEuclidean(G^i); 

12 for vertex w G T.getVerticesInCircle(CE) do 

13 add {v,w) to -E; 


/* Prop. */ 


14 return G\ 


that supports efficient range queries. We adapt the network generation algo¬ 
rithm to use this quadtree and prove subsequently that it achieves subquadratic 
generation time. 


3.1 Generation Algorithm 


Transformation from hyperbolic geometry. Neigh¬ 
bors of a query point u = {4>h,'rh) lie in a hyper¬ 
bolic circle around u with radius R. This circle, 
which we denote as H, corresponds to a Euclidean 
circle E in the Poincare disk. The center Ec and 
radius rad^; of E are in general different from u 
and R. All points on the boundary of E in the 
Poincare disk are also on the boundary of H and 
thus have hyperbolic distance R from u. The points 
directly below and above u are straightforward to 
construct by keeping the angular coordinate fixed 
and choosing the radial coordinates to match the 
hyperbolic distance: and re^), with 

rei,re2 e [0,1), and <iisi-H{Ec,{4>h,re)) 

follows (for details see Appendix [A|): 



Fig. 2: Polar quadtree 


= R for re e {rei,re2}. It 


Proposition!. E^ is at (4>h,-^^) and rad^ is y ~ '^ab+t ’ 

a = cosh(i?) — 1 and b = {1 — r^). 
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Algorithm. The generation of G = (V, E) with n vertices and average degree k 
is shown in Algorithm ^ As in previous efforts, vertex positions are generated 
randomly (lines and We then map these positions into the Poincare disk 
(line and, as a new feature, store them in a polar quadtree (line[^. For 
each vertex v the hyperbolic circle defining the neighborhood is mapped into 


where the neighborhood of v consists of exactly the vertices in the light blue 
Euclidean circle. Edges are then created by executing a Euclidean range query 
with the resulting circle in the polar quadtree (lines [T2p^ . The functions used by 
Algorithm[^are discussed in Appendixj^ We use the same probability distribution 
for the node positions and add an edge (u, v) exactly iff the hyperbolic distance 
between u and v is less than R. This leads to the following proposition: 

Proposition 2. Algorithm^ generates random hyperbolic graphs as defined in 
Section \2. 1[ 


the Poincare disk according to Proposition^ (lines lOpT I - also see Figure lb 


Data structure. As mentioned above, our central data structure is a polar quadtree 
on the Poincare disk. While Euclidean quadtrees are common m , we are not 
aware of previous adaptations to hyperbolic space. A node in the quadtree is 
defined as a tuple (min^, max^, min^., max^) with min^ < max^ and min^. < max^. 
It is responsible for a point p = {fp^rp) S Ili(O) iff (min,^ < fp < rnax^) and 
(miuj. < rp < maxj.). Figure[^shows a section of a polar quadtree, where quadtree 
nodes are marked by dotted red lines. When a leaf cell is full, it is split into four 
children. Splitting in the angular direction is straightforward as the angle range 
is halved: mid,^ := ^ por the radial direction, we choose the splitting 

radius to result in an equal division of probability mass: 


imdru '■= acosh 


cosh(Q! maxr?^) + cosh(Q: min^-H) 


( 3 ) 


(Note that Eq. ^ uses radial coordinates in the native representation, which 
are converted back to coordinates in the Poincare disk.) This leads to two lemmas 
useful for establishing the time complexity of the main quadtree operations: 


Lemma 1. Let Djj be a hyperbolic disk of radius i?, a G K, p a point in D/j 
which is distributed according to Section \2.1\ and T be a polar quadtree on D/j. 
Let C be a quadtree cell at depth i. Then, the probability that p is in C is 4“®. 


Lemma 2. Let R and D/j be as in Lemma^ Let T be a polar quadtree on 
containing n points distributed according to Section 2.1 
large, height(T) € O(logn) whp. 


Then, for n sufficiently 


3.2 Time Complexity 

The time complexity of the generator is in turn determined by the operations of 
the polar quadtree. 


6 










Quadtree Insertion. For the amortized analysis, we consider each element’s initial 
and final position during the insertion of n elements. Let h(T) be the final height 
of quadtree T, let h{i) be the final level of element i and let t{i) be the level of 
i when it was inserted. During insertion of element i, t{i) quadtree nodes are 
visited until the correct leaf for insertion is found, the cost for this is linear in 
t{i). When a leaf cell is full, it splits into four children and each element in the 
leaf is moved down one level. Over the course of inserting all n elements, element 
i thus moves h{i) — t(i) times due to leaf splits. To reach its final position at 
level h{i), element i accrues cost of 0{t{i) + h{i) — t{i)) = 0{h{i)) C 0{h{T)), 
which is O(logn) whp due to Lemma The amortized time complexity for a 
node insertion is then: r(Insertion) S O(logn) whp. 

Quadtree Range Query. Neighbors of a vertex u are the vertices within a Euclidean 
circle constructed according to Proposition Let Af{u) be this neighborhood set 
in the final graph, thus deg(u) := |A/’(u)|. We denote leaf cells that do not have 
non-leaf siblings as bottom leaf cells. A visualization can be found in Figure 

Lemma 3. Let T and n be as in Lemma A range query on T returning a 
point set A will examine at most 0{^/rL -|- |A|) bottom leaf cells whp. 

Due to Lemma the number of examined bottom leaf cells for a range query 
around u is in 0{y/n + deg(u)) whp. The query algorithm traverses T from the 
root downward. For each bottom leaf cell b, 0{h{T)) inner nodes and non-bottom 
leaf cells are examined on the path from the root to b. Due to Lemma h{T) is 
in O(logn) whp. The time complexity to gather the neighborhood of a vertex u 
with degree deg(u) is thus: r(RQ(u)) S O ((i/n -I- deg(u)) • logn) whp. 

Graph Generation. To generate a graph G from n points, the n positions need 
to be generated and inserted into the quadtree. The time complexity of this 
is n ■ O(logn) = O(nlogn) whp. In the next step, neighbors for all points are 
extracted. This has a complexity of 

T(Edges) = ^ O ((-y/n + deg(u)) • logn) = O -b logwhp. (4) 

V 

This dominates the quadtree operations and thus total running time. We conclude: 


Theorem 1. Generating random hyperbolic graphs can be done in + 

m) logn) time whp for sufficiently large n, i. e. with probability > 1 — 1/n. 

4 Experimental Evaluation 

We first discuss several structural properties of networks and use them to analyze 
random hyperbolic graphs generated with different parameters. Comparisons to 
real-world networks and existing generators follow, as well as a comparison of 
the running time to a previous implementation [3j. 
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(a) Clustering coefficient (b) Degree assortativity (c) Max core number 


Fig. 3: Properties of graphs generated with n = 10®, 7 S [2.2, 7] and k S [4, 256]. 
Values are averaged over 10 runs. 


4.1 Network Properties 


We consider several graph properties characteristic of complex networks. The 
degree distribution of many complex networks follows a power law. The clustering 
coefficient is the fraction of closed triangles to triads (paths of length 2 ) and 
measures how likely two vertices with a common neighbor are to be connected. 
Degree assortativity describes whether vertices have neighbors of similar degree. 
A value near 1 signifies subgraphs with equal degree, a value of -1 star-like 
structures. Many real networks have multiple connected components, yet one 
large component is usually dominant. k-Cores are a generalization of components 
and result from iteratively peeling away vertices of degree k and assigning to each 
vertex the core number of the innermost core it is contained in. The diameter 
is the longest shortest path in the graph, which is often surprisingly small in 
complex networks. Complex networks also often exhibit a community structure, 
i. e. dense subgraphs with sparse connections between them. The fit of a given 
vertex partition to this structure can be quantified by modularity m- 

The plots in Figures and (appendix) show the effect of average degree k 
and degree distribution exponent 7 on properties of the generated network. We 
choose a range of [4, 256] for k and [2.2, 7] for 7 , as these cover the ranges of real 
networks in Table Since the geometric model inherently promotes the formation 
of triangles, the clustering coefficient (Fig. 3a) is between 0.6 and 0.9, significantly 


higher than in Erdos-Renyi graphs and thus more realistic for some applications. 
The maximum core number corresponds closely to the density (Fig. 3c).The 
graph is almost alw ays connected for fc > 32 and almost always dis connected 
for fc < 16 (Fig. 6 a). The diameter of the largest component (Fig. 6 b) is highest 
for graphs with fc = 32 and a high 7 , as then the component is large and few 
high-degree hubs exist to keep the diameter small. We use a modularity-driven 
algorithm |20j to analyze the community structure, the size of communities grows 
with the average degree (Fig. | 6 c[ ). Dense graphs with few communities have a 
relatively low modularity (Fig.| 6 e[). Figure [ffi| shows the likelihood ratio of a power 
law fit to an exponential fit of the degree distribution, confirming the power-law 
property. Finally, Figures to [^compare graphs with 10^ vertices generated by 
























our implementation and the implementation of [3]. The differences between the 
implementations are within the range of random fluctuations. 

4.2 Comparison with Real-world Networks 

We judge the realism of generated graphs by comparing them to a diverse set of 
real complex networks (Table [^: PGPgiantcompo describes the largest connected 
component in the PGP web of trust, caidaRouterLevel and as-Skitter repre¬ 
sent internet topology, while citationCiteseer and coPapersDBLP are scientific 
collaboration networks, soc-LiveJournal and fb-Texas84 are social networks 
and uk-2002 and wiki_link_en result from web hyperlinks. Random hyperbolic 
graphs with matching k and 7 share some, but not all properties with the real 
networks. The clustering coefficient tends to be too high, with the exception of 
coPapersDBLP. The diameter is right when matching the facebook graph, but 
higher by a factor of « 100 otherwise, since the geometric approach produces 
fewer long-range edges. Adding 0.5% of edges randomly to a network with average 
degree 10 reduces the diameter by about half while keeping the other properties 
comparable or unchanged, as seen in Figures |10| and [TTI in the appendix. Just 
as in the real networks, the degree assortativity of generated graphs varies from 
slightly negative to strongly positive. Generated dense subgraphs tend to be 
a tenth as large as communities typically found in real networks of the same 
density, and are not independent of total graph size. Similar to real networks, 
random hyperbolic graphs mostly admit high-modularity partitions. Since the 
most unrealistic property - the diameter - can be corrected with the addition of 
random edges, we consider random hyperbolic graphs to be reasonably realistic. 

4.3 Comparison with Existing Generators 

In our comparison with some existing generators, we consider realism, flexibility 
and performance. Typical properties of networks generated by different generators 
can be found in Table The Barabasi-Albert model [5] implements a preferential 
attachment process to model the growth of real complex networks. The probability 
that a new vertex will be attached to an existing vertex v is proportional to 
v's degree, which results in a power-law degree distribution. The distribution’s 
exponent is fixed at 3, which is roughly in the range of real-world networks. 
However, the degree assortativity is negative and the clustering coefficient low. 
The running time is in 0{n^), rendering the creation of massive networks infeasible. 
The Dorogovtsev-Mendes model is designed to model network growth with a 
fixed average degree. It is very fast in theory (0(n)) and practice, but at the 
expense of flexibility. Clustering coefficient, degree assortativity and power law 
exponent of generated graphs are roughly similar to those of real-world networks. 
The Recursive Matrix (R-MAT) model [TD] was proposed to recreate properties 
of complex networks including a power-law degree distribution, the small-world 
property and self-similarity. The R-MAT generator recursively subdivides the 
initially empty adjacency matrix into quadrants and drops edges into it according 
to given probabilities. It has 0(mlogn) asymptotic complexity and is fast in 
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Table 1; Measured properties of some generative models. Parameter ranges are 
n = 10®, fc G [2^, 2®], 7 G [2.2,7] for random hyperbolic graphs, 0.5% additional 
long-range edges for random hyperbolic graphs with long-range edges, n = 
10®, nO G [0,10®), fc G [0,10^) for the Barabasi-Albert generator, PGPgiantcompo, 
caidaRouterLevel, citationCiteseer and coPapersDBLP for Chung-Lu and BTER 
and scale = 16, eF = 10, a G [0.4,1.0), h G [0, a), c G [0, a) for R-MAT. 


Name 

par am 

m 

cc 

deg. ass. 

power law 

7 

diam. 

RHG 

fc,7 

(fc/2) • n 

0.75-0.9 

-0.05-0.7 

yes 

« 7 

3-16k 

RHG-LR 

fc,7 

(fc/2) • n 

0.75-0.9 

-0.05-0.4 

yes 

« 7 

3-30 

Barabasi-Albert 

k, nO 

k ■ n 

0-0.68 

< 0 

if nO < 0.3n 

« 3 

< 30 

BTER 

dd,ccd 

«Edd/2 

matched 

^matched 

possible 

ssmatched 

varies 

Chung-Lu 

seq 

« E seq/2 

< 10"^ 

< 10"^ 

possible 

varies 

8-12 

Dor.-Men. 


2 • n 

0.7 

0.02-0.05 

yes 

5-6 

15-40 

R-MAT 

a, b, c,eF 

eF • n 

0-1 

0-0.6 

yes 

0-10 

0-18 


practice. At least the GraphSOO benchmark parameters[ 6 ] lead to an insignificant 
community structure and clustering coefficients, as no incentive to close triangles 
exists. Given a degree sequence seg, the Chung-Lu (CL) model [T] adds edges (u, v) 
with a probability of p{u,v) = i recreating seq in expectation. The 

model can be conceived as a weighted version of the well-known Erdds-Renyi (ER) 
model and has similar capabilities as the R-MAT model [19]. Implementations 
exist with 0{n -\-m) time complexity [T3|. It succeeds in matching the degree 
distributions of the first four graphs in Table [^ but in all results both clustering 
and degree assortativity are near zero and the diameter too small. 

BTER [12] is a two-stage structure-driven model. It uses the standard ER 
model to form relatively dense subgraphs, thus forming distinct communities. Af¬ 
terwards, the GL model is used to add edges, matching the desired degree distribu¬ 
tion in expectation [IS]. This is done in 0{n-\-mlogd^a,x), where dmax is the maxi¬ 
mum vertex degree. We test BTER with the PGPgicintcompo, caidaRouterLevel, 
citationCiteseer and coPapersDBLP networks. The degree distributions and 
clustering coefficients are matched with a deviation of « 5%. Generated communi¬ 
ties have a size of 5-45 on average, which is smaller than typical real communities 
and those of random hyperbolic graphs. 

As indicated by Table [^ random hyperbolic graphs can match a degree 
distribution exponent, have stronger clustering than the Ghung-Lu and R-MAT 
generator and are more scalable and flexible than the Barabasi-Albert generator. 
Diameter (without additional random edges) and number of connected compo¬ 
nents are less realistic than those produced by BTER, but community structure 
is closer to typical real communities. 

4.4 Performance Measurements 

Figure 1 ^ shows the parallel running times for networks with 10 ^- 10 ^ vertices 
and up to 1.2 • 10® edges. Measurements were made on a server with 256 GB 
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• 

n = 10"^, our impl. 

• 

n = 10^, our impl. 

• 

n = 10®, our impl. 

• 

n = 10^, our impl. 


n — 10^, impl. of [3] 

♦ 

n = 10®, impl. of [3] 


n = lO'^, theoretical fit 


n = 10®, theoretical fit 


n = 10®, theoretical fit 


n = 10^, theoretical fit 



Fig. 4: Comparison of running times to generate networks with 10^-10^ vertices. 
Circles represent running times of our implementation, diamonds the running 
times of the implementation of [3]. Our running times are fitted with the equation 
T[n,m) = ((3.8 • + 1.14 • + 1.38 • 10“®to) logn) seconds. 


RAM and 2x8 Intel Xeon E5-2680 cores at 2.7 GHz. We achieve a throughput 
of up to 13 million edges/s. Even at only 10^ vertices, our implementation is 
two orders of magnitude faster than the implementation of [3] for the same 
graphs. (Note that their implementation supports a more general model.) Due to 
the smaller asymptotic complexity of + m) logn), this gap grows with 

increasing graph sizes. This complexity we prove in Section is supported by 
the measurements, as illustrated by the lines for the theoretical fit. 

5 Conclusions 

In this work we have provided the first generator of random hyperbolic graphs 
with subquadratic running time. Our parallel generator scales to large graphs 
that have many properties also found in real-world complex networks. The main 
algorithmic improvement stems from a polar quadtree, which we have adapted 
to hyperbolic space and which can thus be of independent interest. 

The incremental quadtree construction admits a dynamic model with vertex 
movement, which deserves a more thorough treatment than would have been 
possible here given the space constraints. It is thus part of future work. 

Acknowledgements. We thank F. Meyer auf der Heide for helpful discussions. 
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Appendix 


A Derivation of Proposition 


When given a hyperbolic circle with center {(j)h,rh) and radius rad?t as in Sec¬ 
tion the radial coordinates re of points on the corresponding Euclidean circle 
can be derived with several transformations from the definition of the hyperbolic 
distance: 

vadh = acosh (l + (^) 

2(r 

^ cosh(rad/,) - 1 = \ (6) 

{cosh{radh) - 1){1 - rj) = C^) 


To keep the notation short, we define a := cosh(rad?i) — 1 and 6 := (1 — 
Since radh > 0 and Vh G [0,1), both a and b are greater than 0. It follows: 


(cosh(rad/i) - 1)(1 - r^) = 


2(r2 - 2r/ire r^) 




^ a — a ■ = 


a = r 


2 ‘^{rl-2rhre + rl) 



^a = r^{a+ -) + 


4^0 = r2(a+^)+re^ + ^-a 
b b b 


^0 = ri 


2 r -2 

h 


-drt 

'■b{a+l) ■ b{a+l) (a+f) 


rl)- 

( 8 ) 

(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 


Solving this quadratic equation, we obtain: 


. 2r,, // 2rh V 2r^ - 

ab + 2 \j\ab + 2j ab + 2 


(14) 


Since {tph, rej and ((/f^, Ve^) are different points on the border of E, the center 
Ec needs to be on the perpendicular bisector. Its radial coordinate te^ is thus 
(rei +re^)/2 = To determine the angular coordinate, we need the following 

lemma: 


Lemma 4. Let H be a hyperbolic circle with center at {(phjrh) and radius rad?i. 
The center Ec of the corresponding Euclidean circle E is on the ray from {4>h-,fh) 
to the origin. 
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Proof. Let p be a point in H, meaning dist>^(p, {(ph,rh)) < rad^,. Let p' be the 
mirror image of p under reflection on the ray going through and p. 

The point {4>h,fh) is on the ray and unchanged under reflection: [4>h^fh) = 
{f’hif’h)'■ Since reflection on the equator is an isometry in the Poincare disk 
model and preserves distance, we have dist-uip', {4>h,'rh)) = dist'H(p, [4‘h,'rh)') = 
dist'H(p, {4’h,rh)) < rad/i and p' £ H. The Euclidean circle E is then symmetric 
with respect to the ray and its center must lie on it. □ 

The radius of the circle is then derived from the distance of the center 
to [4>h,ref) and {(fh,re 2 ), which is ~ ■ With both radial and 

angular coordinates of Ec fixed, Proposition follows. 

B Methods Used in Algorithm 

B.l getTargetRadius 

For given values of n, a and i?, the expected average degree k is given by [TSj 
Eq. (22)] and the notation ^ = {a/C,)/{a/C, — 1/2): 

k = len + e-"* („| (^ (i ] _ _ 1)1 + _ 2)) - l)) 

(15) 

As mentioned in Section the value of C, can be fixed while retaining all degrees 
of freedom in the model and we thus assume C = 1 - We then use binary search 
with fixed n, a and desired k to find an R that gives us a close approximation of 
the desired average degree k. 


B.2 mapToPoincare 

In the native representation[T3], the radial coordinate of a point p-^ = r^) 

is set to the hyperbolic distance to the origin: 

r-H = distnipn, (0,0)) 

A mapping g : Ti,^ —>■ Lli(O) from the native representation to the Poincare 
disc model needs to preserve the hyperbolic distance to the origin across models. 
By using the definition of the Poincare metric, we can derive its radial coordinate 
Te in the Poincare disc model: 





I cosh(r^) — 1 
cosh(r^) + 1 


( 16 ) 
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This mapping then gives the correct hyperboiic distance with the Poincare 
metric (Eq. 


dist?^( 5 (( 0 H,r«)), ( 0 , 0 )) = acosh 


1 + 2 


ll5((-^HT«))-(0,0)|p 


(l-||g((0„,r«))|P)(l-||(O,O)|P) 


= acosh 


= acosh 


= acosh 


= acosh 


= acosh 


= acosh 


1 + 2 

( 

1 + 2 

1.2 


{l-\\g((^n,rn)Wm 


'H, 


cosh(r-H)-l ) 112 

cosh(r7^) + l J 11 


( 1 - 


cosh(r.^)-l \ 112', 
y cosh(r-K) + l J ' 


f /cosh{ru)- 


^ y cosh(r,H) + l J 

1 / cosh{r-H)-l\ 

1 y cosh{r-5 

' cosh(r-H)— 1 ) 

K)+iy 

) 

,^cosh(rK) + l ) 

( cosh(r.„)- 



{cosh{r-H) - 1) 


cosh(r^) + 1 — (cosh(r^) — 1 ) 
(cosh(r«) - 1 ) 


(17) 

(18) 

(19) 


( 20 ) 


( 21 ) 

( 22 ) 

(23) 


= acosh ((cosh(r^)) = 


(24) 

(25) 


B.3 transformCircleToEuclidean 

The circle is constructed according to Proposition [l] 


C Proof of Lemma [T] 

To prove Lemma we use Eq. 0 and an auxiliary lemma. 

Lemma 5. Let p he a point in and C a quadtree cell delimited hy min^, 
maxr, min^ and max^. The probability of p being in C is given by the following 
equation: 


Pr(p e C) 


max,^ — min 0 cosh(a max^) — cosh(Q; min^) 
27r cosh(ai?) — 1 


(26) 
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Proof. Let g{{4>, r)) be the probability density for a point p at {fi, r). In Section 2.1 


we used a uniform distribution over [ 0 , 27r) for the angles and defined /(r) as 
the probability density for radial coordinate r. Since the two parameters are 
independent, we write: 


g: [0,27r) x [0,1] 

1 \ 1 sinh(Q;r) 

"(('^’")) = 2^-/(") = 2^-“cosh(ai?)-l 


(27) 

(28) 


With C delimited by min^, max^, min^ and max,^ and g being the product 
of two independent functions, we write: 


^maXr- ^max^ 

Pr(p GC)= / 

J miiir- 'J min^ 


max. /-max.#, sinh(ar) , , , 

• OL -——r- -dodr 

2tt cosh(ai?) — 1 


Constant factors can be moved out of the integral: 


P>-(peC) = 4 


1 


^max.^ ^max^ 


27r cosh(ai?) — 1 
The integrand is independent of (f: 

max0 — min^ 


'mm.^ 


Pr(p e C) = 
Finally, we get: 

Pr(p g C) = 


27r 


a sinh(Q!r)d0(ir 


a sinh(ar)dr 


cosh(ai?) — 1 

max,^ — min^ cosh(a max.) — cosh(a min.) 


27r 


cosh(ai?) — 1 


□ 


We proceed by proving Lemma a by induction. 

Proof. Start of induction (i = 0): At level 0, only the root cell exists and covers 
the whole disk. Since C = D/j, Pr(p g C) = 1 = 4“°. 

Inductive step (i —)■ * + 1): Let Ci be a node at level i. Ci is delimited by the 
radial boundaries min. and max., as well as the angular boundaries min^ and 
max 0 . It has four children at level i + I, separated by mid. and mid^. Let SW 
be the south west child of Ci. With Lemmathe probability of p g SW is: 


Pr(p g SW) = 


mid^ — min^ 


cosh(Q;mid.) — cosh(amin.) 
cosh(Q;i?) — 1 


(29) 


The angular range is halved (mid^ := mid. is selected according 

to Eq. 0: 


mid. := acosh 


/ cosh (a max.) + cosh (a min.) 

V 2 


/a 
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This results in a probability of 


max^ + min^ 
- ^2 - - ~ 

27r 


1 / 1 / coshfa maXr-l+coshfcK miiir-) \ / \ i / • \ 

cosn(a • acosh I- ^^^^ 1 /a) — cosn(amin^) 


cosh{aR) — 1 


(30) 


„ . f cosh(amax^)+cosh(amin^) \s. 

max^ + min 0 —2 min^ cosn(^acosn 2 " ' J > 

47 r cosh(Q;i?) — 1 

. cosh(Q; maXr-)+cosh(a mirir) i / • \ 

max0 — min0 -^-^2-^- - ~ cosh(a mm^j 

47 r cosh(ai?) — 1 

max 0 — min 0 cosh(a max^) + cosh(a minr) — 2 cosh(Q! minr) 
47 r 2 (cosh(aii) — 1 ) 

1 max 0 — min^ cosh(a max^) — cosh(a minr) 

4 27r cosh(ai?) — 1 

= ipr(pGC 0 


cosh (a min^) 

(31) 

(32) 

(33) 

(34) 

(35) 


As per the induction hypothesis, Pr(p G Ci) is 4“* and Pr(p G SW) is thus 
I • 4“® = Due to symmetry when selecting mid,^, the same holds for the 

south east child of Ci. Together, they contain half of the probability mass of Ct. 
Again due to symmetry, the same proof then holds for the northern children as 
well. □ 


D Proof of Lemma 

We say “with high probability” when referring to a probability of at least 1 — 1/n 
(for sufficiently large n). While previous results exist for the height and cost of 
two-dimensional quadtrees these quadtrees differ from our polar hyperbolic 
approach in important properties and the results are not easily transferable. For 
example, we adjust the size of our quadtree cells to result in an equal division of 
probability mass when taking the hyperbolic geometry into account, see Lemma 
We thus make use of a lemma from the theory of balls into bins instead: 

Lemma 6 ( |15j L When n balls are thrown independently and uniformly at ran¬ 
dom into n bins, the probability that the maximum load is more than 3 Inn/ In Inn 
is at most 1 /n for n sufficiently large. 

Proof (of Lemma^. In a complete quadtree, 4*^ cells exist at height k. For 
analysis purposes only, we construct such a complete but initially empty quadtree 
of height k = [log 4 (n)], which has at least n leaf cells. As seen in Lemma a 
given point has an equal chance to land in each leaf cell. Hence, we can apply 
Lemma with each leaf cell being a bin and a point being a ball. (The fact 
that we can have more than n leaf cells only helps in reducing the average load.) 
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From this we can conclude that, for n sufficiently large, no leaf cell of the current 
tree contains more than 3 Inn/In Inn points with high probability (whp). Even 
if we had to construct a subtree below a current leaf I to store points whose 
number exceeds the capacity of I, the height of this subtree cannot exceed the 
number of points in the corresponding area, which is at most 3 Inn/In Inn whp. 
Consequently, the total quadtree height does not exceed O(logn) whp. 

Let T' be the quadtree as constructed in the previous paragraph, starting 
with a complete quadtree of height k and splitting leaves when their capacity is 
exceeded. Let T be the quadtree created in our algorithm, starting with a root 
node, inserting points and also splitting leaves when necessary, growing the tree 
downward. 

Since both trees grow downward as necessary to accomodate all points, but 
T does not start with a complete quadtree of height fc, the set of quadtree nodes 
in T is a subset of the quadtree nodes in T'. Consequently, the height of T is 
bounded by O(logn) whp as well. □ 

E Proof of Lemma [3] 

Proof. As done previously, we denote leaf cells that do not have non-leaf siblings 
as bottom leaf cells, see Figurefor an example. The following proof is done for a 
leaf capacity of one. Since a larger leaf capacity does not increase the tree height 
and adds only a constant factor to the cost of examining a leaf, this choice does 
not result in a loss of generality. 

Let L be the set of bottom leaf cells containing a vertex in A and let Q be 
the set of bottom leaf cells examined by the range query. Since the contents of 
leaf cells are disjoint, \L\ < |A| holds. The set Q\L consists of leaf cells which 
are examined by the range query but yield no points in A. These are empty leaf 
cells within the query circle as well as cells cut by the circle boundary. 

Empty leaf cells occur when a previous leaf cell is split since its capacity is 
exceeded by at least one point. Therefore an empty leaf cell a in the interior of 
the query circle only occurs when at least two points happened to be allocated 
to its parent cell b. A split caused by two points creates four leaf cells, therefore 
there are at most twice as many empty bottom leaf cells as points. 

The number of cells cut by the boundary can be bounded with a geometric 
argument. On level k = |"log 4 n], at most cells exist, defined by at most 2^ 
angular and 2^ radial divisions. When following the circumference of a query 
circle, each newly cut leaf cell requires the crossing of an angular or radial division. 
Each radial and angular coordinate occurs at most twice on the circle boundary, 
thus each division can be crossed at most twice. With two types of divisions, the 
circle boundary crosses at most 2 • 2 • 2* = 4 • cells on level k. Since the 

value of 4 • 2 r*°S 4 "l is smaller than 4 • this yields < 8 • ^/rl cut cells. 

In a balanced tree, all cells on level k are leaf cells and the bound calculated 
above is an upper bound for |(3\L|. For the general case of an unbalanced tree, 
we use auxiliary Lemma[^ which bounds to 0{^/n) the number of bottom leaf 
cells descendant from cells cut in level k. □ 
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Fig. 5: Visualization of bottom leaf cells in a quadtree. Bottom leaf cells are 
marked in green, non-bottom leaf cells in red and interior cells in black. 


Lemma 7. Let T, n, R, Djj and a be as in Lemma^ Let k := |"log 4 n~\ and let 
C he a set of [c ■ ^/n\ quadtree cells at level k, for c > 1. The total number of 
bottom leaf cells among the descendants of cells in C is then bounded by 4c • y/n 
whp. 


Proof. New leaf cells are only created if a point is inserted in an already full 
cell. We argue similarly to Lemma that the descendants contain at most twice 
as many empty bottom leaf cells as points. The number of points in the cells 
of C is a random variable, which we denote by X. Since each point position is 
drawn independently and is equally likely to land in each cell at a given level 
(Lemma [^, X follows a binomial distribution: 

X~B („, (36) 

For ease of calculation, we define a slightly different binomial distribution Y ~ 
B{n, hV"). Since n < 4^ and c • -^/n > [c • ^/n\ , the tail bounds calculated for Y 
also hold for X. 

Let H (be the relative entropy (also known as Kullback-Leibler 

divergence) of the two Bernoulli distributions B{ and B{‘^^^). We can 

then use a tail bound from to gain an upper bound for the probability that 
more than 2c points are in the cells of C: 


Pr(V > 2c • -/n) < exp 





(37) 


For consistency with our previous definition of “with high probability”, we 
need to show that Pr(y > 2c^pn) <lln for n sufhciently large. To do this, we 
interpret Pr(T > 2c • ^/rl)/{\/n) as an inhnite sequence and observe its behavior 
for n —>■ oo. Let a„ := Pr(y > 2c • ^/rl)/{l/n) = n ■ Pr(y > 2c • ^/n) and 


bn ■= n • exp 



From Eq. (371 we know that a„ < b„ 


19 










Using the definition of relative entropy, we iterate over the two cases (point 
within C, point not in C) for both Bernoulli distributions and get: 


bn = n ■ exp ( —nH 


2c\/n c^/n 


2 c 


= n ■ exp —n —j= In 2 + 1- -j= In 


m 


2 c 


m 


1 - 


2c^/n ' 


1 - 


2 c 

/n 


= n ■ exp I —n—^ In 2 ) • exp ( —n ( 1- — ) In 

/n — c 


2 c 

fn 


Cy/n I I 
n / / 

\Jn — 2c 


= n ■ exp (—2c-\/nln2) • exp ({n — 2cy/n^ In ^ 

\ — 2c J 


y/n - c 
22cV" Vv^-2c 


n—2c\/n 


4c-v/n 


1 + 


y/n — 2c 


'i—2cy/n 


(38) 

(39) 

(40) 

(41) 

(42) 

(43) 


(While hn is undefined for n S {c^,4c^}, we only consider sufficiently large n 
from the outset.) 

1 

We apply a variant of the root test and consider the limit lim„_>oo(6„) vs for 
an auxiliary result: 


lim 

n—¥oo 




= lim n vs . _ 
n—>-(x> 4 ^ 





^/n — 2c 


(44) 

(45) 

(46) 


From e/4 < 0.7, c > 1 and the limit definition, it follows that almost all 
1 

elements in (6„) vs are smaller than 0.7 and thus almost all elements in are 
smaller than 0.7^^. Thus lim„_>oo bn < lim„_>oo 0.7^ = 0. Due to Eq. (37), we 
know that an is smaller than bn for large n, and therefore that the number of 
points in C is smaller than 2c • y/n with probability at most f/n for n sufficiently 
large. Again with high probability, this limits the number of non-leaf cells in C to 
c • y/n and thus the number of bottom leaf cells to 4c • y/n, proving the claim. □ 
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F Further Graph Property Parameter Studies 


- 7 = 2.2 
- 7 = 4.6 
7 = 7.n 


103 

■i lO-* 
I 10® 
& 10^ 
10 ' 


#- 7 = 2.2 
•- 7 = 4.6 
• 7 = 7.0 


(a) Vertices in largest com- (b) Diameter of largest (c) Avg. community size 
ponent component 


- fc = 4 

- fc = 32 
• fc = 256 


*- 7 = 2.2 
•- 7 = 4.6 
» 7 = 7.0 


7 = 2.2 
7 = 4.6 
-«• 7 = 7.0 


(d) Measured vs. desired 7 (e) Modularity of commu- (f) Loglikelihood of power- 
nities. law degree distribution 


Fig. 6: Further parameter studies, omitted from Figure]^ due to space constraints. 
Values are averaged over 10 runs, except for the diameter. 
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G Properties of Some Real Networks 


Table 2; Properties of various real networks. The columns show the number of 
vertices and edges, the clustering coefficient, the maximum core number, the 
log likelihood of a power-law degree distribution, the exponent of an optimal 
power-law fit, the degree assort at ivity, the diameter, average size of communities 
and modularity of the community structure. 



n 

m 

cc 

max core 

pi 11 

7 

deg.ass. 

diameter 

comm, size 

mod. 

PGPgiantcompo 

lOK 

24K 

0.44 

31 

2.04 

4.41 

0.23 

24 

101 

0.88 

fb-TexasS4 

36K 

1.6M 

0.19 

81 

1.54 

4.8 

0 

7-8 

1894 

0.38 

caidaRouterLevel 

192K 

609K 

0.19 

32 

6.73 

3.46 

0.02 

26-30 

365 

0.85 

citationCiteseer 

268K 

IM 

0.21 

15 

9.6 

3.0 

-0.05 

36-40 

1861 

0.80 

coPapersDBLP 

540K 

15M 

0.81 

336 

4.04 

5.95 

0.50 

23-24 

2842 

0.84 

as-Skitter 

1.7M 

IIM 

0.3 

111 

20.3 

2.35 

-0.08 

31-40 

1349 

0.83 

soc-Live Journal 

4.8M 

43M 

0.36 

373 

6.94 

3.34 

0.02 

19-24 

632 

0.75 

uk-2002 

18.5M 

261M 

0.69 

943 

331 

2.45 

-0.02 

45-48 

441 

0.98 

wiki_link_en 

27M 

547M 

0.10 

122 

26 

3.41 

-0.05 

- 

21.6 

0.67 


H Comparison with Previous Implementation[3J 

Both implementations sample random graphs, making a direct comparison of 
generated graphs difficult. In its output files, the implementation of [3] provides 
the hyperbolic coordinates of the generated points. Yet, since the distance 
threshold R is computed non-deterministically with a Monte Carlo process and 
not written to the log file, we do not have all necessary information to recreate 
the graphs exactly. The Figures mill and show properties of the generated 
graphs instead, averaged over 10 runs. Plots showing graphs created with the 
implementation of [3] are on the left, plots created with our implementation are 
on the right. Some random fluctuations are visible, but for almost all properties 
the averages of our implementation are very similar to the implementation of 
[3]. The measured values of 7 for thin graphs and various target 7 s differs from 
the previous implementation, but the fluctuation within the measurements of 
each implementations are sufficiently strong that it leads us to assume some 
measurement noise. The differences between the implementations are smaller 
than the variations within one implementation. 
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Clustering Coefficient 
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Fig. 7: Comparison of clustering coefficients, size of largest component and diam¬ 
eter of largest components for the implementation of [3] (left) and our implemen¬ 
tation (right). Values are averaged over 10 runs. 
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Degree Assortativity 
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Measured vs Desired 7 
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Fig. 8 : Comparison of degree assortativity, degeneracy and measured vs desired 7 
for the implementation of [3] (left) and our implementation (right). Values are 
averaged over 10 runs. 
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Loglikelihood of Powerlaw Fit 


Loglikelihood of Powerlaw Fit 



Fig. 9: Comparison of likelihood of a power-law fit to the degree distribution 
for the implementation of [3] (left) and our implementation (right). Values are 
averaged over 10 runs. 


I Effect of Additional Random Edges 
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Fig. 10: Effect of random long-range edges on diameter ranges. Baseline graphs 
are generated with an average degree of 10, values are averaged over 5 runs. Black 
circles correspond to PGPgiantcompo, caidaRouterLevel, citationCiteseer and 
as-Skitter from Table [2l which have a comparable density. 
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(a) Clustering coefficient. (b) Loglikelihood of power-law degree dis¬ 

tribution. 



Fig. 11: Side effects of adding random edges to generated graphs. Baseline graphs 
are generated with an average degree of 10, values are averaged over 5 runs. The 
clustering coefficient changes by less than 0.03, the two lines of the likelihood of 
a power-law degree distribution are nearly identical. 
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